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ABSTRACT 


In order to develop a computer program able to model the 
frictional heating of metals in high pressure oxygen or 
nitrogen a number of additions have been made to the 
frictional heating model originally developed for tests in low 
pressure helium. These additions include: l) a physical 

property package for the gases to account for departures from 
the ideal gas state, 2) two methods for spatial discretization 
(finite differences with quadratic interpolation or orthogonal 
collocation on finite elements) which substantially reduce the 
computer time required to solve the transient heat balance, 3) 
more efficient programs for the integration of the ordinary 
differential equations resulting from the discretization of 
the partial differential equations and 4) two methods for 
determining the best-fit parameters via minimization of the 
mean square error (either a direct search multivariable 
simplex method or a modified Levenburg-Marquardt algorithm) . 
The resulting computer program has been shown to be accurate, 
efficient and robust for determining the heat flux or friction 
coefficient vs. time at the interface of the stationary and 
rotating samples. 
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INTRODUCTION 


The frictional heating apparatus at WSTF is used to determine 
the flammability of various metals undergoing frictional 
heating in oxygen-rich environments. A detailed description of 
this apparatus has been given elsewhere 1 as well as a list of 
the previous mathematical models for the frictional heating. 

As a part of the 1991 NASA/ASEE Summer Faculty Fellowship 
Program a computer model 2 of the frictional heating apparatus 
was developed for frictional heating experiments in low 
pressure helium (approximately 75 psia) . The model that was 
developed was successful in accurately fitting the 
experimental data for the helium via a best— fit of the 
friction coefficient or heat production at the interface of 
the rubbing samples. 

The success of the model for the simplest case (i.e., a low 
pressure non— reactive gas) led to the hope that, similar 
results could be achieved in much more difficult situations 
(i.e., a high-pressure reactive gas). The goal of this project 
was therefore to expand the low pressure model to include high 
pressure oxygen and nitrogen systems (at least up to 1,000 
psia) . 


MODEL DEVELOPMENT 


The model used in this study is the same transient one- 
dimensional equation for heat conduction as in the original 
model. 2 For each annular sample (i = 1/2) we have the 

following energy balance 




8 k dT -Q -Q ~Q 
TZ* 1 ?* Vv Uc Ur 


( 1 ) 


where = the density of sample i, C pi — the heat capacity of 
sample i, T = temperature in the sample, t = time, x = axial 
position in the sample, k, = the thermal conductivity of sample 
i, q v = the heat loss per unit rod volume due to convection, 
Q = the heat loss per unit rod volume due to conduction and 
Q r = the heat loss per unit rod volume due to radiation. 


The initial condition is: 

T = T a for -L 2 <> x < L t (2) 

where L, - the length of the each sample (i = 1 is the 
stationary sample and i = 2 is the rotating sample) . 
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The boundary conditions are: 
At x = 0: 


~^i + ^2 " Qf ant ^ T 0 , - T 0 _ ( 3 ) 

where Q F is the total energy flux at the interface via 
friction and/or reaction. 

At x = L t (i = 1/ 2) : 


T - T w - T a (4) 

Energy Flux at the Interface: 

In a nonreactive environment the energy flux Q F is given by 
the equation 

C f (t) = P(t)v(t) fi(t) in Btu/in 2 sec. (5) 

The applied pressure (P) and velocity (v) terms are measured 
quantities but the coefficient of friction (n) must be 
determined via a best-fit of the experimental data. For cases 
where both friction and reaction are important it is possible 
to simply use the heat flux Q F as the parameter to be fit. 
Three correlations are possible for fitting p or Q p vs. time: 

1. The first model simply states that the parameter to be 

fit, p(t) (e.g., the heat flux Q F or friction coefficient 

H) , will be held constant over a specified time interval 
At,, i.e. 

p(t) = p, for t, <> t <> t,+At, (6) 

2. The second model uses an Nth-degree polynomial to 

approximate the parameter p(t) via N+l coefficients, 
i.e. , 

p(t) = a 0 + a, t + a 2 t 2 + + a N t N (7) 

3. The third model uses parameters similar to Model 1 but 
also uses quadratic interpolation to approximate the 
value of the parameter p(t) as a continuous function of 
time for all times t, £ t £ t,+At,. The parameter values 
used for interpolation are p, through p 4+2 . 
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Gas-Phase Heat Balance; 


For the gas phase the following energy balance was used 
(assuming the gas phase to be well-mixed) : 

m g C V ' 9 ^ = 2nR 0 h v (T - T g )dx - - T v ) (8) 

where m = the mass of gas in the system, C„ ^ = the constant 
volume heat capacity of the gas, h v — convective heat transfer 
coefficient, T = the rod temperature (varies with position) , 
T = the temperature of the gas, T w — the chamber wall 
temperature, R 0 = the outer rod radius and A g = the area for 
heat transfer from the gas to the surrounding chamber . 

Due to the high pressures (above 1,000 psia) used for the 
oxygen and nitrogen systems, the assumption of ideal gas was 
inappropriate. This was especially true for physical 
properties such as heat capacity, thermal conductivity and 
viscosity. It was therefore decided to use the Lee— Kesler 
modification of the Benedict— Webb— Rubin eguation of state 3 to 
approximate the departure from the ideal gas properties. (The 
derivation of the departure equations was carried out by a 
graduate student (K. Bhattacharrya) under NASA Grant NAG9-557 
and will be discussed in more detail in his thesis. It is 
important to note that the equations for the departure 
function for gas heat capacity given by both reference 3 and 
the original paper of Lee and Kesler are incorrect.) 

External Heat Transfer From Rods: 

The following equations are used for the heat transfer losses 
from the rods: 

Convection : The heat loss from the rods due to convection can 
be described by the eguation 

Q v = SjlJT - T t ) (9) 

where S v = the surface area per unit rod volume available 
for convective heat transfer (= 2 RJ (R„-R?) ) and R, = the 
inner radius of the annulus. The convective heat transfer 
coefficient is assumed to be given by a Nusselt number 
correlation as: 


Nu = K(Ta/F g ) a Pr b 


( 10 ) 
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where K, a and b are constants given by one of three 
methods: 1) the correlation of Becker and Kaye 4 , 2) the 
correlation of Eisenberg et al. 5 or 3) a best-fit of the 
data. The modified Taylor number is defined as 

u?r d 3 

(Ta/F ) ■ — " ( 11 ) 

vF t 

where o> = the angular velocity of the rotating cylinder, 
“ {R c + R 0 )/2, d - R c - R 0 , R c = the radius of the 
chamber, v = the kinematic viscosity of the gas and F = 
a shape factor given by Becker and Kaye 4 . 


For heat conduction from the rods to the gas the Nusselt 
number is defined as 


Nu m 


2h v R 0 


( 12 ) 


where k g = the thermal conductivity of the gas. For heat 
conduction from the gas to the chamber the Nusselt number 
is defined as 


Nu * 



(13) 


Conduction: For the heat transfer from the cylinders to the 
holders the heat loss due to conduction is approximated 
by a •'pseudo" heat transfer coefficient h c 


Qc - *cA(r - TJ 


S cf t 

Ax 


(T - T w ) 


(14) 


where * the surface area per unit volume of rod i 
available for conductive heat transfer and Ax ■= the 
distance separating the cylinder and holder. 

Radiation : The heat loss due to radiative heat transfer is 
calculated via 


Q r = S r /JF H (T 4 - T 4 ) (15) 

where all temperatures are absolute (°R), a *= Stefan's 
constant * = the surface area per unit volume of rod 

i available for radiative heat transfer and Fy = the 
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shape factor for radiative heat transfer from body i to 
body j . 


NUMERICAL SOLUTION 


The addition of the physical property package for the gases at 
high pressures resulted in a significant increase in the 
computation time for the model. It therefore became apparent 
that a more efficient method for the spatial discretization as 
well as the time integration would be required to keep the 
computation times to realistic values. 

For the time integration the original subroutine LSODE was 
replaced by LSODES (which consists of the same integration 
methods but was designed to more efficiently handle the sparse 
Jacobian matrices which result from reducing partial 
differential equations to a series of ordinary differential 
equations via spatial discretization) . LSODES was written by 
A. C. Hindmarsh and A.H. Sherman. The Runge-Kutta-Fehlberg 
program RKF45 of H.A. Watts and L.F. Shampine was also 
included in the program because of its increased robustness. 
(Both LSODES and RKF45 were obtained through the NETLIB 
distribution system.) 

For the spatial discretization two different methods were 
incorporated. One method is the same finite difference method 
as in the original model but, in order to determine the 
temperatures at the exact axial thermocouple locations, 
quadratic interpolation was used when these positions were 
located between grid points. This allowed a significant 
reduction in the required number of grid points (from 18 to 5) 
to achieve the same accuracy with respect to measured vs . 
calculated temperatures. 

The second method of spatial discretization was to use 
orthogonal collocation on finite elements. 6 This method 
consists of replacing the original partial differential 
equation with the ordinary differential equation 

At x = x n : 

' = A £ AqTj - Q v (T h ) - Q c {T h ) - Q r (T„) 

at h 2 Pi pi 

where x„ = the nth collocation point in the element, h 
size of the element, A, v « the elements of the (N+2) x 


(16) 

= the 
(N+2) 
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matrix approximating the first derivative via orthogonal 
collocation 6 and N = the number of interior collocation points 
in each element. Finite elements were used to simplify the 
integration of the convective heat transfer contribution in 
the gas energy balance (which was accomplished via 
quadrature) . 

Two elements were used in each rod with the boundary of each 
element being the point where the sample holders ended. The 
solution at the collocation point at the boundary between the 
two elements in each sample was handled in the usual way by 
setting the axial gradients equal to each other. If we number 
the four elements sequentially with element 1 being the 
section of the stationary sample in the holder we get the 
following equation for the boundary of the elements in the 
stationary sample 


h 



N°>* 2 


E 


Vi" - 


_1 


N**2 

E 


where x tJ = the axial position for the location of the boundary 
between the elements in the stationary sample, T® = the 
temperature in element j next at the collocation point i in 
each element. For the rotating sample we simply replace the 
superscript (1) with (4) and (2) with (3) and x t , with x er , the 
axial position for the element boundary in the rotating 
sample . 


The boundary condition at the interface was handled in a 
similar manner 


A-l 


A 


N+2 

E a m t “ * Of 


e,r *•! 


To make the program more robust two parameter fitting 
subroutines have been added. One method is a direct search 
multivariable simplex method developed by Nelder and Mead. 7 
The other subroutine is LMDIF which uses a modification of the 
Levenburg-Marquardt algorithm and was written by B.S. Grabow, 
K.E. Hillstrom and J.J. More. (LMDIF was obtained through the 
NETLIB distribution system. ) 


Results 


The computer program (now called FHP and installed on the 
computer system at WSTF) was tested on five high pressure 
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systems , all with oxygen at approximately 1,000 psia, and gave 
the results shown in Table 1. The error in Table 1 refers to 
the root mean square error between the measured and calculated 
values for both thermocouples in the stationary sample, 0.05 
and 0.20 inches from the interface. Tests 830-198 and 830-199 
used AISI 304 stainless steel and tests 830-305 through 830- 
307 used AISI-316 stainless steel. Physical properties of 
these materials vs. temperature were obtained from Touloukin. 


TABLE 1. RESULTS FOR HIGH-PRESSURE OXYGEN TESTS 


Test 

Average Error 
( °F) 

830-198 

±26.4 

830-199 

±63.6 

830-305 

±14.3 

830-306 

±18.8 

830-307 

±19.4 


The temperatures in the samples varied from initial 
temperatures of approximately 70 °F to final measured 
temperatures of from 1,000 to 1,400 °F. 

Figure 1 gives an example of the results of using Model 3 to 
approximate the energy flux at the interface for test 830-307. 
Thermocouple TC-702 is located 0.05 inches from the interface 
of the rotating and stationary samples and TC-703 is 0.20 
inches away. Also shown is the predicted interface 
temperature. For the AISI 304 tests (830-198 and 830-199) the 
calculated temperature drop between TC-702 and TC-703 was 
significantly smaller than the measured value. This was 
probably due to inaccurate data for the physical properties of 
the metals. 

Figure 2 shows a comparison of the heat flux (in Btu/in 2 sec) 
at the interface (via Model 3) vs. time (in seconds) for the 
three repeated tests 830-305, 830-306 and 830-307. 

Conclusions 

The program FHP, developed as a result of the support over the 
two past s umm ers by the NASA/ASEE Summer Faculty Fellowship 
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Program, is successful in accurately determining the heat flux 
and/or friction coefficients necessary to produce the observed 
temperatures in the WSTF frictional heating apparatus. The 
program has enough alternatives with respect to spatial 
discretization, integration and minimization to make it very 
robust and sufficiently accurate compared to possible 
inaccuracies in the experiments. The weak points in the 
program are: 

1. The quantitative effects of convection on the heat 
loss from the samples. At present a simple 
correlation from studies on annuli with the inner 
cylinder rotating is used. 

2. The physical properties of some of the metals vs. 
temperature may be inaccurate. 

3. The assumption that the holders are at a constant 
temperature throughout the simulation. 

4. The one-dimensional nature of the model. 

At present, only the last problem is actively being 
investigated. As a part of the research under NASA grant NAG9- 
557 a computer model is being developed in order to determine 
the possibility of significant radial gradients through the 
samples for heat fluxes at the interface which vary radially 
as well as with time. Once this is known a model can then be 
developed to include angular variations in the heat flux to 
determine whether there is a significant difference in 
temperatures at the same axial position but different angular 
positions . 
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Results for Model 3 
[Experiment 830-307] 



Figure 1. Comparison of measured and calculated 
temperatures for thermocouples located 
0.05 inches (TC-702) and 0.20 inches (TC- 
703) from the interface. (T in °F and t 
in sec.) 
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Heat Flux for Model 3 



t 

830-305 

830-306 

830-307 


Figure 2. Heat flux at the interface ( Q F in 
Btu/in 2 sec.) vs. time (in sec.) for three 
different tests using Model 3 . 
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